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Abstract. We study simultaneous price drops of real stocks and show that for high drop thresholds they fol- 
low a power-law distribution. To reproduce these collective downturns, we propose a minimal self-organized 
model of cascade spreading based on a probabilistic response of the system elements to stress conditions. 
This model is solvable using the theory of branching processes and the mean-field approximation. For a 
wide range of parameters, the system is in a critical state and displays a power-law cascade-size distribution 
similar to the empirically observed one. We further generalize the model to reproduce volatility clustering 
and other observed properties of real stocks. 
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1 Introduction 



Cascade spreading is an important emergent property of 
various complex systems. Real life examples of cascades 
are numerous and range from infrastructure failures and 
epidemics to traffic jams and cultural fads JTJ[2] - Theoret- 
ical models of cascades usually assume that agents can 
be in one of two states (healthy or failed) and an agent's 
failure puts some stress on its neighbors which may con- 
sequently fail too. See [3] for a recent survey of this field 
offering a novel unifying view. 

In this paper we focus on cascades in economic sys- 
tems which can be identified with stock prices suddenly 
dropping in a major market crash [3] or with companies 
going bankrupt simultaneously and leading to global re- 
cession [5] . Theoretical models of such cascades are based 
on shortage and bankruptcy propagation in production 
networks [B] , default propagation in credit networks p3[5] , 
interaction of firms through one monopolistic bank [S] or 
in a complex credit network economy |10j . and herding 
behavior of traders [TTl[T2"] . While these models help us to 
understand cascade processes in economic systems, they 
are mostly too involved to allow for analytical solutions — 
their study hence relies on numeric simulations and agent- 
based modeling [P5] . 

A simpler point of view on cascade phenomena is of- 
fered by the concept of self-organized criticality (SOC) 
which has had a deep impact on the science of complex- 
ity. First introduced more than twenty years ago to ex- 
plain the ubiquitous 1// noise [2], it caused a blossoming 
of toy models, computer simulations, and real life exper- 
iments |15j . The analytical techniques employed include 
scaling arguments [16], mean-field theories [17], branch- 
ing processes Q2], renormalization methods (191120) . and 
rigorous algebraical techniques [21] ■ 



SOC is a mechanism which explains the emergence of 
complex behavior in many diverse real world systems [221 
|2"5] . The generic behavior of SOC models is: (a) they evolve 
so that they always stay close to the critical point, (b) 
long periods of robustness and moderate activity are inter- 
rupted by sudden breakdowns. This qualitatively resem- 
bles "stock markets which expand and grow on relatively 
long time scales but contract in stock-market crashes on 
relatively short time scales" [TS] and "stock crashes caused 
by the slow buildup of long-range correlation leading to a 
global cooperative behavior of the market eventually end- 
ing into a collapse in a short time interval" [4] . This simi- 
larity provides the main motivation for the present study. 

We begin our work with an empirical investigation 
of simultaneous price drops of real stocks and show that 
the size distribution of observed events is broad (for high 
drop thresholds it follows a power-law distribution) . This 
observation suggests that simultaneous stock downturns 
are a collective phenomenon. We propose a simple dy- 
namical model which for a wide range of parameters self- 
organizes into a critical state. Unlike most SOC models, 
our model assumes a probabilistic response mechanism 
where a node has only a certain probability of reacting to 
the current stress conditions. The basic idea behind mod- 
eling simultaneous stock downturns with cascades is that 
decline of a single stock may provoke investors' reactions 
which consequently may cause other stocks to decline and 
a "cascade" to spread. The key premise is that while failed 
nodes become significantly more resistant in the next time 
step, healthy nodes become slightly less resistant. This 
close parallel with the slow growth/fast decay picture de- 
scribed above is further supported by our analysis of em- 
pirical data which shows that majority of stocks behave 
in this way. While there are certainly many other effects 
contributing to the dynamics of market crashes (external 
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shocks, for example), we show that failure propagation 
alone can reproduce some of the observed patterns. 

The minimal model proposed here has the advantage 
of being simple, not relying on fine-tuning of parameters, 
analytically solvable in some cases, and easily generaliz- 
able to more complicated settings. We analyze it using the 
formalism of branching processes, the mean-field approxi- 
mation and, for complex topologies of nodes' interactions, 
using numerical simulations. Obtained cascade-size distri- 
butions exhibit a close similarity to our empirical obser- 
vations. Introduction of memory within the model allows 
us to reproduce other empirically observed features, such 
as volatility clustering, though at the cost of analytical 
tractability. We conclude our study with a discussion of 
further model's generalizations and possible areas of ap- 
plication. 



2 Empirical data 

Here we investigate co-occurring price movements of real 
stocks. Adopting the vocabulary of cascade models, we say 
that a stock fails when the relative loss of its price over 
a given time interval At exceeds a certain threshold H. 
Denoting the price of stock i at time t as Pi(t), its failure 
occurs when [pi(t) — Pi{t + At)]/pi(t) > H. The number 
of stocks failing at time t, np(t), is a direct analog of the 
cascade size in a model of cascade spreading. As the input 
data we use daily closing prices (hence At = 1 day) of 500 
stocks from the standard U.S. index S&P 500 (this data 
is freely available at, for example, finance.yahoo.com). 
To achieve a fixed system size, we consider only those 
332 companies which are in the stock market since the 
beginning of 1992 and use their prices during the 18-years 
long period ending in May 2010 for our analysis. 

The empirical distribution of failure sizes is shown in 
Fig. [1] for H = 0% and H = 10%. We see that for the 
large value of H (which is in line with the notion of stock 
failures), the observed size distribution has a power-law 
shape. Using the methodology described in [51], we ob- 
tained the power-law exponent 2.19 ± 0.05 with the lower 
bound for the power-law behavior n m j n = 3. The corre- 
sponding p- value (obtained using the standard Kolmogo- 
rov-Smirnov statistic) is 0.92 which confirms that the data 
is consistent with the hypothesis of a power-law distribu- 
tion. Similar results are obtained also for other threshold 
values so long as H > 8%. When H < 8%, the resulting 
size distributions are broad but probably not power-law. 
Finally, when H — 0% (i.e., any price drop is interpreted 
as a failure), the size distribution is roughly symmetric 
around the value corresponding to one half of the system 
size (see Fig.[T]). In the following analysis of empirical data 
we use the threshold H = 10%. 

The power-law shape itself suggests that the observed 
simultaneous stock downturns are rather a collective phe- 
nomenon than independent events. This hypothesis is fur- 
ther supported by the average correlation of simultane- 
ously failing stocks, 0.35 (again including only events with 
at least three simultaneously failing stocks), which is sig- 
nificantly higher than the overall average stock correla- 




Fig. 1. The empirical failure size distribution observed with 
real stock prices (daily closing prices of 332 companies from 
January 1992 until May 2010) for threshold relative drops 
H = 0% and H — 10%. The straight line corresponds to the 
exponent 2.19 obtained by statistical analysis of the data. 



tion, 0.25. Another sign of a strong connection among si- 
multaneously failing stocks comes from their division to 
ten different industrial sectors according to the GICS clas- 
sification. The effective number of sectors participating in 
a cascade is defined as 



10 

i=l 



(1) 



where is the relative share of sector i in the cascade and 
2i=i r i = 1- By averaging this quantity over all cascades 
of a given size S, we obtain e(S). This number can be com- 
pared with the effective number of sectors corresponding 
to selecting failed stocks at random, e'(S). The analysis of 
stock prices shows that for any S > 3, e(S) is significantly 
smaller than e'(S) which implies that simultaneous stock 
failures preferentially affect strongly connected stocks in 
one sector or in a small number of sectors. 

Now we turn our attention to time correlations of fail- 
ures. The autocorrelation of the number of failing stocks 
with the time lag one day, C(nF(t),np(t + 1)) 0.15, is 
comparable with the autocorrelation of absolute returns, 
C(\r(t)\, \r(t + 1)|) « 0.25 (the latter result agrees with 
previous studies [25,26 ). The positive autocorrelation val- 
ues are signs of volatility clustering which is commonly ob- 
served in financial data [27]. (Loosely speaking, volatility 
clustering means that large changes tend to be followed 
by large changes and small changes tend to be followed by 
small changes, as first noted by Mandelbrot [28].) 

We further estimate conditional failure probabilities 
for individual stocks. For example, P(F\N) denotes fail- 
ure probability of a stock given that this stock didn't fail 
in the previous time step (other three quantities, P(N\F), 
P(F\F), and P(N\N), follow the same logic). When the 
results are averaged over all stocks, we obtain P(F\F) = 
0.039 which is much higher than the overall failure proba- 
bility P(F) — 0.003 — this is another sign of volatility clus- 
tering in our data. On the level of individual stocks, how- 
ever, 62% of all stocks with at least three failures strongly 
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satisfy the inequality P(N\F) > P(N) which is equiva- 
lent to P(F\F) < P(F) (because P(F\F) + P{N\F) = 1). 
(By strong satisfying we mean that the difference of the 
two probabilities is greater than the sum of their uncer- 
tainties.) We see that despite volatility clustering in the 
data, most stocks are more "resistant" to failures after 
they have just undergone one. For the remaining stocks, 
probabilities P(N\F) and P{N) either differ less than the 
sum of values' uncertainties (for 14% of stocks) or even 
strongly satisfy the opposite inequality P(N\F) < P(N), 
with corresponding values of P(F\F) often as high as 0.30 
(24% of stocks). 

To summarize, after a failure (a major price drop), 
most stocks become more resistant to another failure — 
this observation will serve as a basis for the mathemati- 
cal model presented in the following section. At the same 
time, there is a fraction of stocks which are prone to con- 
secutive failures — this particular feature will be discussed 
in detail in Section 2] 



3 Basic model and its mean-field solution 

In this section we present a basic model which is amenable 
to analytical treatment and qualitatively reproduces some 
of the features observed in empirical data. In its original 
formulation, this model is particularly suitable for stocks 
that, as discussed in the previous section, after a failure 
become more robust. A generalization of the model aiming 
at reproducing other observed features (volatility cluster- 
ing, for example) is presented in Section 2] 

Consider a system of N nodes where node i (i — 
1, . . . , N) has only two possible states: failed (i € J 7 ) and 
healthy (i F). With each node i we further associate 
fragility fi 6 [0, 1] which measures how this node reacts 
to failures of its neighbors (the higher the fragility, the 
more likely is the node to follow a neighbor's failure). The 
dynamics of the model is governed by the following sim- 
ple rules, (i) In each time step, the first failed node ("trig- 
ger" ) is chosen at random and may induce failures of other 
nodes, (ii) If a neighbor of node i fails, node i follows it 
with probability /$ and resists with probability 1 — fi. (If 
several neighbors of node i fail simultaneously, in order 
to stay healthy, node i has to resist each individual fail- 
ure.) The cascade of failures propagates until all remaining 
nodes resist the damage, (iii) At the end of the time step, 
fragilities of all nodes are updated according to 



Mt + i) = { 



(i + 



i e T 



(2) 



where < /3 <C 1 and A 6 (0, 1) are parameters of 
the model (in effect, failed nodes become less fragile and 
healthy nodes become slightly more fragile in the next 
time step). All values fi(t + 1) > 1 are truncated to 1 
(this may occur when j3 is large). After this update is 
finished, all nodes are again marked as healthy, the cur- 
rent time step ends and a new one begins with point (i). 
Note that unlike some other models of cascade spread- 
ing, failed nodes are not removed from the system in our 



case. If a long enough equilibration period is applied be- 
fore measuring the system behavior, the initial fragility 
values fi(0) are of little importance (see Section [3~51 for a 
detailed discussion). Unless stated otherwise, we set them 
randomly in the range (0, 1) in our simulations. 

According to the rules above, when n neighbors of node 
i fail, node i resists with the probability (1 — fi) n and fails 
with the complementary probability 



P F (f u n) = l-{l-f i ) n . 



(3) 



This response to failures is "path-independent" in some 
sense: the probability that a node resists n failures of its 
neighbors, (1 — fi) n , is the same as the probability of re- 
sisting two consequent waves of failures of x and n — x 
neighboring nodes, (1 — fi) x {l — fa) n ~ x - 

We simplify the system by assuming that interactions 
of all nodes are equally strong (the general case will be 
studied in Section I3~4")l . This renders the notion of "node's 
neighbors" superfluous because every failure affects all re- 
maining healthy nodes in the system. Now assume that 
after the initial failed node is chosen, ri\ nodes respond to 
this failure and fail too. Each of the remaining N — hq — ni 
nodes (here no = 1 is the initial number of failed nodes) 
then has some ni-depcndcnt failure probability which re- 
sults in 712 new failures, and so on, until in iteration m, 
n m = is achieved. The cascade size is then defined as 
the total number of failures, S = uq + • • • + n m , and node 
fragilities are consequently updated according to Eq. 
Since in one turn nodes can only fail once, cascade sizes 
are limited by the system size and S < N. 

The dynamics of the system, based on failure propaga- 
tion and fragility updating, is fully contained in the three 
above-described rules. In the following paragraphs we shall 
study when these rules drive the system to a critical state 
and what is the distribution of cascade sizes P(S). 

3.1 Failure probability 

Let Pp be the average failure probability of a given node 
in one time step (or, equivalently, the average fraction 
of failed nodes in one time step). Assuming that t cq is 
some sufficiently long equilibration time (we use t cq = 
10 4 for all our simulations), later fragility values averaged 
over realizations, (fi), do not evolve anymore. All nodes 
interact equally strongly, hence (fi) is independent of i 
and it can be replaced with (/). Since in a large number 
of time steps T each node undergoes PfT failures and 
(1 — Pf)T non- failures, Eq. @ implies 



(f(t cq + T)) = (f(t cq ))X p - T (l + P) 



(1-Pp)T 



(4) 



Using the equilibrium condition (f{t cq + T)) = (/(t eq )), 
we can solve this equation with respect to Pp to get 



ln(l + /3) 
In 



(5) 



When j3 <C 1, this can be approximated with Pf(/3, A) s» 
— /3/lnA (Fig. [5] compares these results with numerical 
simulations). 
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Fig. 2. Average failure probability: Pf given by Eq. (O (solid 
black line), Pf ~ — /3/lnA (dashed blue line) and numerical 
results (symbols, averaged over 10 6 time steps) for N = 10 3 , 
A = 0.1. The vertical dotted line indicates Pq given by Eq. ([6]). 



As already mentioned, when interactions of all nodes 
are equal, (fi) is independent of i. If wc further neglect 
fluctuations of fi, then all nodes have identical fragility 
(/) . This is a mean-field-like approximation which replaces 
the exact cascade spreading with cascade spreading in a 
homogeneous averaged medium. Since the number of di- 
rect descendants now follows a simple binomial distribu- 
tion with mean N(f), we can use elementary results of 
branching process theory to express the average cascade 
size (the total progeny) as (S) = 1/(1 — N{f)). Further, 
using (S) — NPp(f3, A) we obtain the average fragility 



(/> = 



1 - 



In [(l + /3)/A] 
AHn(l + (3) 



(7) 



Since /3 > and A < 1, (/) is always less than 1/N. Com- 
parison with numerical simulations (not shown) confirms 
that Eq. © is valid only for /3 < 1. 



A node may fail because it is selected as the first failed 
node (with probability 1/N) or due to failure propagation 
(with probability Pp); Pf thus can be written as Pf = 
1 /N+Pp. Since the value of Pf depends solely on /3 and A, 
Pp = Pp — 1/N may be negative for a small system which 
is, of course, impossible in practice. This situation occurs 
when for given A, N, the value of (3 is smaller than a certain 
threshold (3q and hence it docs not suffice to compensate 
for the fragility decay due to A. Eq. (HJ) then has only 
the trivial solution (/) = and hence Pf(J3,X) = 1/N 
(failures do not spread). When j3 is small, the approximate 

form of Pf can be used to solve this equation with respect where p a is defined using 
to and we get 

A) « -— (6) 



3.3 Cascade size distribution 

The theory of branching processes is well studied [31] and 
can be easily applied to our model. According to a theo- 
rem from [32] . if the generating function for the number 
of direct descendants d is 7r(x), the total progeny of the 
resulting branching process Y has the distribution 



N 

which agrees with numerical simulations (see the vertical 
line in Fig. [2]). Note that if the number of initial failed 
nodes is assumed to grow with the system size as wN 
(w <C 1), we get fio ~ —win A which is independent of N. 

When model parameters are set to extreme values (for 
example, N = 10 3 , /3 = 10 3 , A = 10~ 3 ), the system ex- 
hibits unusual modes of behavior where active turns (with 
nearly all nodes failed) alternate with calm turns (with 
nearly all nodes healthy) . While Eq. ((5|) holds also in such 
conditions, our further analysis focuses on j3 <C 1 which 
renders more realistic behavior. 



3.2 Average fragility 

When nfi <C 1 , Pf given by Eq. Q can be approximated 
as Pp(/i,n) ~ nfi which can be interpreted as indepen- 
dence of stress inflicted by n individual failed nodes. This 
further means that each failed node has its failing descen- 
dants independently of other failed nodes and hence one 
can use the theory of branching processes to describe 
the cascade spreading. Note that by use of this theory we 
implicitly assume that the system size is infinite. For a dis- 
cussion of the finite-size effects on the size of an epidemic 
outbreak see [50] , 



P(Y\n ) = ^ P ?l no 



[n(x)] b =p { b) + pf ) x + 



(8) 



(9) 



and no is the number of ancestors (in our case, the number 
of initial failed nodes) . Since d obeys a binomial distribu- 
tion, its generating function is tt(x) = (1 — (/) + (f)x) N 
and we get 



P(S\(3,X) = 



1 ( NS 



S\S 



_J</> 



s-i 



[i-(f)f s ~ s+1 (10) 



where we used no = 1 and (/) is given by Eq. ([7]). Note 
that the resulting probability is positive for S > N which 
contradicts the model assumptions (each node fails at 
most once in a given turn). This is a direct consequence 
of using the theory of branching processes which assumes 
that the system size is infinite. This problem is of little im- 
portance for small values of (3 when the obtained values 
of P(S) are negligible for S > N. 

When 1 <C S <C N, Eq. (TIT))) can be approximated with 



P(S\0,\) 



(N(f)) 



S-i S(l-N(f)) 



/2nS 3 / 2 



(11) 



According to Eq. ([7]), limAr^oo N(f) = 1 for any given /3, A 
and hence in the limit of large system size is P(S\(3, A) ~ 
g -3 / 2 which corresponds to the classical critical branching 
process. For a finite system, the smaller the value of /3, the 
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larger the value of 1 — N(f). Consequently, the power-law 
scaling holds only for S <C f3N (this agrees with Fig. [3] 
where for f3 = 10~ 3 , the power-law behavior disappears 
at S ~ 10). On the other hand, the range of (3 and A for 
which the system self-organizes to a critical state is wide 
and we can say that this is an SOC system. 

A comparison of the obtained analytical results with 
numerical simulations is shown in Fig. [3] The agreement 
is good for small values of /3 (/3 < 0.01) and the initial 
slope of the distributions (before the finite-size effects be- 
come apparent) is close to —3/2. Results obtained with 
(3 = 0.001 confirm that when (3 is small enough, P(S) 
decays faster than as a power law. When (3 is large, true 
P(S) deviates from the analytical prediction and exhibits 
a secondary maximum at a large size value — this effect is 
well visible in Fig. [3] for /3 ~ 0.1. This maximum, formally 
simply a super-critical phase of the model, resembles so- 
called meaningful outliers discussed in '.'>■'> . To estimate 
the value of (3 at which the secondary maximum appears 
and Eq. (fTU|) ceases to hold, we take the average number 
of failures computed both from Eq. ([TO]) and from Eq. ([5]). 
By comparing the two results we obtain 

N 

NP F ((3,X) = J2SP(S\(3,X). (12) 
s=i 

When (3 is small, both sides of this equation depend on 
j3 and the equality can hold. However, Eq. (jllj) shows 
that when f3 is sufficiently large, the size distribution is 
approximately power-law and it is independent of f3. As 
we increase (3 further, the power-law distribution does not 
suffice to provide enough failures and for Eq. ([T2"j) to hold, 
an additional contribution must appear on the rights side. 
The value /3i when this happens can be found by substi- 
tuting P(S) ~ S*~ 3 / 2 on the right side and approximating 
the summation with integration. When N is large, we ob- 
tain 

/ 2 \ 1/2 

ft«-(^J In A (13) 

which complements the previously found threshold /3q . For 
N = 10 4 and A = 0.1, we obtain f3\ w 0.02 which agrees 
with our empirical observation (/? < 0.01 for Eq. (fTUf to 
hold) above. 

Finally, by comparing the empirical observations pre- 
sented in Fig. Q] with the obtained analytical results, we 
can conclude that the presented model exhibits qualitative 
agreement with the studied real system. 

3.4 Generalizations 

To test how robust are the obtained results, we consider 
simple generalizations of the proposed model. First of all, 
when the multiplicative fragility update rule Eq. © is 
replaced by an additive one, the behavior of the system 
does not change considerably. The second generalization 
relates to the assumed even influence of a node's failure on 
all the remaining nodes. Denoting the strength of failure 
propagation from node i to node j as Cij , the probability 




Fig. 3. The cascade size distribution: numerical results (color 
lines), analytical results according to Eq. (1101) (dashed lines) 
and the power-law decay with exponent —3/2 (thick solid line) 
for N = 10 4 , A = 0.1, 10 7 time steps, and /3 = 0.001 (red line, 
fastest decay), /3 = 0.01 (green line, medium decay), /3 = 0.1 
(blue line, slowest decay). The analytical solution is not plotted 
for ft = 0.1 because it is very similar to that for /3 = 0.01. 

that node j fails as a result of Vs failure can be generalized 
to Ci.jfj. The probability that node j fails as a result 
of a group T of failed nodes (given by Eq. J3J) before) 
generalizes to the form 

P F {f j ,T) = l-J{{l-C i>j f j ). (14) 

Matrix C encodes the structure of the network of node 
interactions. 

When the elements Cjj are drawn independently from 
a given distribution and the system size is large, the mean- 
field approximation is again appropriate to describe the 
system behavior and the power-law size distribution with 
exponent 3/2 results. Similarly when C contains a block 
structure with inter-block elements drawn from a differ- 
ent distribution than intra-block elements (this mimics the 
sector structure of the stock correlation matrix [21)1131] ). 
the original power-law size distribution remains largely 
unchanged (unless either the block division of C or one of 
the two probabilistic distributions are such that they do 
not allow to use the mean-field approximation) . Analogous 
behavior results from the "random neighbor approxima- 
tion" in which node's neighbors are chosen anew repeat- 
edly (see [33] for this kind of analysis of a different model). 

When all elements Cj j are either zero or one, matrix C 
can be represented by a network and a complex topology 
of node interactions can be introduced by network mod- 
els [33]. We studied two different types of networks: the 
Erdos-Renyi network where Ci_j = 1 with probability p 
and Ci.j = otherwise and the growing Barabasi-Albert 
network where each new node is attached to / old nodes. 
(These two kinds of networks are structurally very distinct 
as the former consists of nodes of approximately identical 
degree and the latter exhibits a power-law degree distribu- 
tion.) Numerical results for both cases are shown in Fig. 2] 
As expected, for the Erdos-Renyi network with p > 1/-/V, 



6 



Stanislao Gualdi, Matiis Medo, Yi-Cheng Zhang: Self-organized model of cascade spreading 



P(S) 10 



P(S) |() s 




Fig. 4. The cascade size distribution on complex networks: 
(a) sparse Erdos-Renyi networks with p = 5 ■ 10 -5 , 10 -4 , 10 -3 ; 
the indicative thick line has slope 1.5, (b) Barabasi- Albert net- 
works with I — 1 and I = 10; the indicative thick line has slope 
1.65). Parameters of the system: N = 10 4 , /3 = 0.005, A = 0.1, 
10 7 time steps. 



the size distribution exponent remains unchanged. When 
p < 1/N, the network consists of small isolated compo- 
nents and hence big cascades cannot occur. The irregular 
size distribution P(S) observed for j3 = 5 • 10~ 5 is due to 
topological properties of the particular network realization 
where the model was simulated (i.e., positions of respec- 
tive ups and downs of the size distribution depend on the 
network realization). These results agree with a previous 
study of the sandpile dynamics [37] (see [35] for an exten- 
sive recent review of critical phenomena in complex net- 
works). By contrast, Barabasi- Albert networks yield cas- 
cade size distributions with significantly higher exponents 
(approximately 1.65) which is probably due to strong in- 
homogeneity of the network. When 1=1, P(S) deviates 
from a power law, probably as a consequence of the scale- 
free network topology (the same shape of the distribution 
is observed for different realizations of the network). 



3.5 Role of the initial fragility values 

While it sounds plausible that due to model's stochastic- 
ity, the initial fragility values have no influence on the 
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Fig. 5. Fragility distributions at different time steps (the initial 
fragility values are set to 1/N, N = 10 4 , f3 = 0.001, A = 0.1). 



equilibrium fragility distribution, the situation is in fact 
more complicated. For example, a simple numerical sim- 
ulation with fi(0) = 1/N for all i shows a case where: 
(i) no stationary fragility distribution arises, (ii) at any 
time step, only a small number of distinct fragility val- 
ues is observed (see Fig. [5]). What causes the discreteness 
of fragility values? Denoting the number of failing and 
healthy time steps of node i as Fi and Hi , respectively, it 
must hold that Fi + Hi = t where t is the current time 
step. This node's fragility now can be written as 



f i (t) = f i (0)(l+/3) t [\/(l+P)Y 



(15) 



When all /i(0) are identical, the possible values of fi(t) 
are discrete at any time step t and the ratio of neighbor- 
ing possible values is (1 + f3)/X. If A is small (as it is in 
our simulations), this ratio is large and hence the number 
of actually observed fragility values is small (because val- 
ues much smaller or greater than the average fragility are 
unlikely). Eq. (fTSj) implies that possible fragility values de- 
pend on t and hence there can be no stationary fragility 
distribution — this is confirmed by Fig. [5] where fragility 
peaks constantly shift to higher values and change their 
relative heights. Interestingly, even this peculiar setting of 
fi(0) does not alter the long-term model's behavior sub- 
stantially and the aggregate quantities (such as the av- 
erage failure probability or the cascade size distribution) 
are similar to those found for randomized initial fragility 
values before. 

Differences between neighboring peaks are A/(l + /3), 
hence the time after which the fragility distribution pat- 
tern repeats can be estimated as In [A/ (1 + /3)] / ln(l + j3). 
Since this is a typical time of fragility evolution, one can 
use it also as an estimate of the initial equilibration time 
T cq . For the smallest value of /? in our simulations (/3 — 
0.005) we obtain T cq ps 4 600 which ex post confirms 
our setting of the equilibration time to 10 4 . Finally, note 
that while the random setting of fi(0) prevents discrete 
fragility values from appearing, some remnants of the ini- 
tial fragility values can be preserved by Eq. (fT5|) . To ob- 
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tain a fragility distribution truly independent of the initial 
values, one has to assume annealed dynamics, i.e. fragility 
updating by randomized values of ft and A. 



4 Generalized model with partial memory 

Fragility updating rules defined by Eq. @ imply that 
nodes become more robust after a failure and hence au- 
tocorrelation of their failures as well as autocorrelation of 
the total number of failures are negative (their magnitudes 
depend on ft and A). As discussed in Section^ this is true 
for majority of stocks but certainly not for all of them. 
To allow for repeatedly failing stocks, we introduce the 
probability a with which a failed node stays failed also in 
the next time step (and consequently acts as an additional 
initial failed node). This probability has the role of par- 
tial memory in the system and, as we shall see, gives rise 
to volatility clustering and other effects observed in real 
financial data. Note that memory or delayed stress propa- 
gation are quite often part of cascade spreading models as 
in, for example, [35]. We assume that fragilities of nodes 
which stay failed due to a are not updated in the given 
time step (when a is small, this assumption has little in- 
fluence on the results). 

Since P(F) <C 1, the probability of a node's repeated 
failure is now P(F\F) w a which, in combination with 
the empirical results presented in Section [21 motivates 
us to set a — 0.04. We further choose ft — 0.01 and 
A = 0.1 which best correspond to the critical regime in 
Fig. [3] Using this setting we numerically obtain condi- 
tional probabilities consistent with those observed in the 
empirical data: P(F\F) = 0.041 (empirical value is 0.039), 
P(F) = 0.004 (empirical value is 0.003) and P{F\N) = 
0.004 (empirical value is 0.003). As long as we stay in the 
critical regime, these values depend on ft and A weakly. 
Presence of volatility clustering is confirmed by signifi- 
cantly positive autocorrelation of the number of failures 
C(nF{t),np(t + 1)) 0.3 (empirical value is 0.15). The 
precise value depends on a and A (and much less on ft) 
but positive autocorrelation naturally arises for a which is 
large enough. By contrast, a = yields P(F\F) ss P(F) 
and C(n F (t),n F (t + 1)) « -0.04. Finally, Fig. H shows 
P(S) for different values of a. We see that for small values 
of a, the size distribution remains power-law with expo- 
nent gradually decreasing as a grows. Due to the addi- 
tional complexity introduced by partial memory, an ana- 
lytical cascade size distribution for this generalized model 
has not been obtained yet. 



5 Discussion 

We studied empirical stock prices and found that large 
simultaneous downturns follow a broad distribution con- 
sistent with a power law with exponent 2.19 ± 0.05. To 
reproduce this behavior, we proposed a minimal stochas- 
tic model of failure propagation. Using a mean field ap- 
proximation and branching process theory we derived the 




Fig. 6. Cascade size distribution for the modified model: nu- 
merical results for ft — 0.01, A = 0.1, N = 10 4 , 10 7 time steps, 
and various values of a. 



general cascade size distribution and determined the range 
of parameters which give rises to the critical regime. To 
reproduce other features observed in financial data, such 
as volatility clustering, partial memory was introduced 
within the basic model. 

While our model implicitly assumes arrival of news to 
the market (they cause the initial nodes to fail and allow 
cascades to be created), we minimize the influence of news 
on the system's behavior by assuming their equal impact 
(in each time step, exactly one initial node is chosen to 
fail). This approach is motivated by the extensive study 
of excess volatility which shows that it is difficult to link 
the observed trading volumes and volatility to the arriving 
information [3D] and even the large crash of 1987 does not 
seem to be triggered by particular news [41) . In reality, 
of course, the impact of news on the market differs from 
one day to another. It could be therefore interesting to 
test how different ways of choosing the initial failed nodes 
influence the model's behavior (for example, the number 
of the initial nodes can be random or network hubs may 
be preferentially chosen to trigger a cascade). 

There is a number of other challenging questions which 
deserve further investigation. Firstly, since the cascade 
sizes corresponding to the secondary maximum in Fig. [3] 
are comparable with system size, this behavior cannot 
be described within the formalism of branching processes 
where an infinite system size is assumed. While we found 
an approximate condition for the appearance of the sec- 
ondary maximum, how to proceed further towards an ex- 
haustive description of the resulting size distribution is 
still an open question. Secondly, it would be interesting to 
find an analytical expression for the size distribution expo- 
nent in scale-free networks where it appears to differ from 
the mean-field value 3/2. Thirdly, generalized model with 
"partial memory" , studied only numerically here, calls for 
analytical approaches. Fourthly, it would be interesting 
to know whether the model can be modified to produce 
power-law size distributions with exponents considerably 
higher than those reported here. One opportunity for such 
a generalization is to assume a dynamic network structure 
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whose evolution depends on nodes' failures, similarly to 
the approach used in [421143] for different models. Alter- 
natively, as a generalization of the current binary model 
where nodes are either healthy or failed, one could define 
a multi or continuous-state model in which the probability 
of following a neighbor's failure depends on the failure's 
magnitude. 

We stress that the probabilistic spreading mechanism 
proposed here is a general one and its use is not limited to 
market crashes or firm bankruptcies. For example, eco- 
nomic exchanges between countries are so intense that 
decline in one country may propagate to a neighboring 
one (take, for example, how growth in many European 
countries depends on spending of German consumers) . On 
a two or three dimensional lattice, a similar mechanism 
might be employed to model earthquakes because, simi- 
larly to the proposed model, a failure at one place of the 
Earth's crust exerts some stress on its neighborhood (the 
number of failed nodes would then represent the earth- 
quake size). In summary, the proposed model, together 
with its generalizations, has proven to be simple yet rich 
in behavior. It poses a variety of new research questions 
and we are looking forward to its future development and 
applications. 
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